%% US Public Debt and Safe Asset Market Power
%% Jason Choi, Rishabh Kirpalani, and Diego Perez
%% Sept 1, 2023

%% Calibrate Monopoly Equilibrium

disp('% Epsilon = 2.2, Lambda = 1') 
P.gamma = 2;
P.lambda = 1;
rspread = 0.0062;
rus = 0.0053;
busus = 0.41;
P.eta = 0.545;

P.beta = 1/(rspread+(1+rus));
P.nu = rspread/(busus^(P.eta-1));
P.omega = P.nu*P.eta*busus^(P.eta-1-P.lambda);

P.alpha = 0.3;
P.delta = 0.1;
P.delta_rw = P.delta;
P.delta_us = P.delta;

P.gdp_ratio = 1.1;
P.NFA_us = -0.05;
P.home_bias = 0.8;

[param_out err]=fsolve(@(x) func_moments(x,P),[0.7,0.9,1,1]);
P.sigma = param_out(1);
P.sigma_star = param_out(2);
P.A = param_out(3);
P.Astar = param_out(4);


disp('% Epsilon = 1.5, Lambda = 1') 
P.gamma = 2;
P.lambda = 1;
rspread = 0.0062;
rus = 0.0053; %% real rate
busus = 0.41;
P.eta = 0.313;
P.N = 1;

P.beta = 1/(rspread+(1+rus));
P.eta_f = rspread/(busus^(P.eta-1));
P.lambda_f = P.eta_f*(1+1/P.N*(P.eta-1))*busus^(P.eta-1-P.lambda);

P.alpha = 0.3;
P.delta = 0.1;
P.delta_rw = P.delta;
P.delta_us = P.delta;

P.gdp_ratio = 1.1;
P.NFA_us = -0.05;
P.home_bias = 0.8;

[param_out err]=fsolve(@(x) func_moments(x,P),[0.7,0.9,1,1]);
P.sigma = param_out(1);
P.sigma_star = param_out(2);
P.A = param_out(3);
P.Astar = param_out(4);


disp('% Epsilon = 3.2, Lambda = 1') 
P.gamma = 2;
P.lambda = 1;
rspread = 0.0062;
rus = 0.0053; %% real rate
busus = 0.41;
P.eta = 0.690;
P.N = 1;

P.beta = 1/(rspread+(1+rus));
P.eta_f = rspread/(busus^(P.eta-1));
P.lambda_f = P.eta_f*(1+1/P.N*(P.eta-1))*busus^(P.eta-1-P.lambda);

P.alpha = 0.3;
P.delta = 0.1;
P.delta_rw = P.delta;
P.delta_us = P.delta;

P.gdp_ratio = 1.1;
P.NFA_us = -0.05;
P.home_bias = 0.8;

[param_out err]=fsolve(@(x) func_moments(x,P),[0.7,0.9,1,1]);
P.sigma = param_out(1);
P.sigma_star = param_out(2);
P.A = param_out(3);
P.Astar = param_out(4);


disp('% Epsilon = 2.2, Lambda = 0') 
P.gamma = 2;
P.lambda = 0;
rspread = 0.0062;
rus = 0.0053; %% real rate
busus = 0.41;
P.eta = 0.545;
P.N = 1;

P.beta = 1/(rspread+(1+rus));
P.eta_f = rspread/(busus^(P.eta-1));
P.lambda_f = P.eta_f*(1+1/P.N*(P.eta-1))*busus^(P.eta-1-P.lambda);

P.alpha = 0.3;
P.delta = 0.1;
P.delta_rw = P.delta;
P.delta_us = P.delta;

P.gdp_ratio = 1.1;
P.NFA_us = -0.05;
P.home_bias = 0.8;

[param_out err]=fsolve(@(x) func_moments(x,P),[0.7,0.9,1,1]);
P.sigma = param_out(1);
P.sigma_star = param_out(2);
P.A = param_out(3);
P.Astar = param_out(4);


disp('% Epsilon = 2.2, Lambda = 2') 
P.gamma = 2;
P.lambda = 2;
rspread = 0.0062;
rus = 0.0053; %% real rate
busus = 0.41;
P.eta = 0.545;
P.N = 1;

P.beta = 1/(rspread+(1+rus));
P.eta_f = rspread/(busus^(P.eta-1));
P.lambda_f = P.eta_f*(1+1/P.N*(P.eta-1))*busus^(P.eta-1-P.lambda);

P.alpha = 0.3;
P.delta = 0.1;
P.delta_rw = P.delta;
P.delta_us = P.delta;

P.gdp_ratio = 1.1;
P.NFA_us = -0.05;
P.home_bias = 0.8;

[param_out err]=fsolve(@(x) func_moments(x,P),[0.7,0.9,1,1]);
P.sigma = param_out(1);
P.sigma_star = param_out(2);
P.A = param_out(3);
P.Astar = param_out(4);
P
